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Abstract 

We derive semiclassical contributions of periodic orbits from a boundary integral equation 
for three-dimensional billiard systems. We use an iterative method that keeps track of the 
composition of the stability matrix and the Maslov index as an orbit is traversed. Results are 
given for isolated periodic orbits and rotationally invariant families of periodic orbits in axially 
symmetric billiard systems. A practical method for determining the stability matrix and the 
Maslov index is described. 
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1 Introduction 



Billiards are popular models for the study of dynamical systems and their quantum counter- 
parts. They are simpler than general systems with a potential and exhibit a large spectrum of 
dynamical behaviour ranging from integrability to chaoticity. Until recently, most investigations 
have concentrated on two-dimensional billiards whose numerical treatment requires much less 
effort than higher-dimensional systems. But two-dimensional systems also have rather special 
properties like the division of the energy surface into separate regions by invariant tori. 

During the last three years there have been several publications on three-dimensional billiard 



systems as models for chaotic systems 



In other fields three- 



II, I, i i S I, I, H, ITT], H 

dimensional billiard systems have served as models much longer, like in nuclear or clusters 
physics. An overview on these applications is given in |T3[. In this context, Balian and Bloch 



investigated in detail semiclassical approximations for the level density in three-dimensional 
billiard systems |1T4|, O, Il6| . Their results for the contributions of periodic orbits are expressed 



in terms of determinants of 2n-dimensional matrices, where n is the number of reflection points 
of a trajectory. In contrast to this, the results of Gutzwiller for smooth three-dimensional 
systems are expressed in terms of the 4 x 4-stability matrix [0, and the relation between both 
results is not explicit. 

For two dimensional billiard systems Harayama and Shudo derived semiclassical contri- 
butions of isolated periodic orbits in terms of the stability matrix starting from a boundary 
integral equation [|18|. Their derivation involved the transformation of n-dimensional matrices, 
where again n is the number of reflections of a trajectory. For higher dimensional systems this 
method would be very elaborate. In addition, the relation between the obtained index in the 
trace formula and the number of conjugate points is not direct [IT^ . 

We avoid this difficulty by applying an iterative method that follows the trajectory from 
reflection to reflection and keeps track of the composition of the stability matrix and the Maslov 
index as the trajectory is traversed. This method has proven to be convenient for the derivation 



of uniform approximations for diffractive orbits in two-dimensional billiard systems [20 



The motivation for the present work is two-fold: firstly, we want to show how the contribu- 
tions of isolated periodic orbits in terms of the 4 x 4-stability matrix can be obtained directly 
from a boundary integral equation without involving large matrices. Secondly, we want to 
present a practical method for the determination of semiclassical contributions of isolated peri- 
odic orbits and families of orbits in systems with axial symmetry, including the determination 
of the stability matrix and the Maslov index. Furthermore, the method is convenient for the 
treatment of more complex contributions of periodic orbits like those of bifurcations that are 
of interest in nuclear physics applications E^ . 



2 The boundary integral equation 

A powerful method for investigating quantum billiards, numerically as well as analytically, 
consists in replacing the Schrodinger equation and boundary conditions by an integral equation. 
This integral equation is defined on the billiard boundary and effectively reduces the dimension 
of the problem by one. For three-dimensional billiards with Dirichlet boundary conditions the 
integral equation has the following form |2^ 

dadn'Go{r,f',E)dn^{r} = ^9^,.^(f') . (1) 
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Here T) denotes the three-dimensional domain of the bilhard system, and dV is the boundary 
surface. For simphcity we assume in the following that the domain V is compact and the 
boundary dV is . The vectors r and r ' are located on the boundary dV, and da is an 
area element at r. Furthermore, n is a unit normal vector to the boundary pointing outwards, 
and c?fi = n ■ V is the component of the gradient in direction of n. Dimensionless units with 
h = 2m = 1 are used throughout, ip is a. solution of the Schrodinger equation with Dirichlet 
boundary conditions, and Gq denotes the free outgoing Green function which is given by 



G,ir,f',E) = - ^y^-'" , (2) 



expjifc 


r — r ' } 


Ait 


r 







where k = yE. Eq. (|I|) is a Fredholm equation of the second kind. It has non-trivial solutions 
only if the determinant A{E) := det(l — Q{E)) vanishes, where Q is the integral operator 
defined on the boundary dV by 

Qu{f') = -2l dadn,Go{f,f',E)u{r). (3) 

JdV 

In recent years there were several publications on the Fredholm determinant A{E) for two- 
dimensional billiards. The integral equation in two dimensions is non-singular if the boundary 



is continuously differentiable, i.e. if it doesn't have corners or cusps (see |2^, gjj). For the 
following discussion we assume that V is compact and simply connected, and that the boundary 
&D is C^, although the following results are valid under less restrictive conditions. By an 
application of Fredholm theory, A(i?) can be represented by an absolutely convergent series 
which, when approximated semiclassically, yields a series in terms of so-called pseudo-orbits 



I^, Considered as a function of k the determinant A(£') is holomorphic in the k- 



plane except for a cut which is due to the cut of the free outgoing Green function in two 
dimensions. It is known that in the upper half plane Im A; > the determinant has only zeros 
which correspond to the eigenvalues of the Schrodinger equation for the interior billiard with 
Dirichlet boundary conditions k = i^E^. For Im/c < there can be further zeros. These 
additional zeros correspond to the resonances of the outside scattering problem with Neumann 
boundary conditions [^, [28[] . 

The spectral staircase for the interior problem can be represented in the following form 

N{E) = NUE) - - hm Im log ^^^j^ , (4) 

where the second term on the right-hand side jumps by one at every eigenvalue of the Schrodin- 
ger equation (or by the multiplicity of an eigenvalue if it is degenerate). Ns^{E) is a smooth 
function given by 

1 A' 

NsmiE) = - / dx \im lm—(x + ie) . (5) 



It has been conjectured in and shown in that TiNsm{E) is equal to the total scattering 
phase ?7n(-E') = | logdet S'n(-E) for the outside scattering problem with Neumann boundary 
conditions. The on-shell scattering matrix for this problem is denoted by Sf^{E). 

The asymptotic expansions of the smooth parts N{E) of the spectral staircase and f]-!^{E)/7r 
of the total scattering phase agree in the first three leading terms and in general differ in the 
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next terms |2^. A relation between the complete asymptotic expansions of the smooth part 
N[E) for the inner billiard problem and of the smooth part f]{E) for the outside scattering 
problem with the same boundary conditions was conjectured and proved for several leading 
terms (13 terms for the Dirichlet case) in It is stated that if the asymptotic expansion of 
N{E) is 

oo 

N{E) ^aE + hE^'^ + c + Y^ c„E^/2-n ^ 



n=l 



then the asymptotic expansion of ri{E) is 

n{E) 



aE-hE^I^ + c-Y,CnE^'^-^ . (7) 



^ 1 



A semiclassical evaluation of the second term on the right-hand side of @ yields contribu- 
tions from interior periodic orbits as well as from periodic orbits that run in the exterior of the 
billiard (if they exist) [jl9|, The contributions from the exterior orbits have to be cancelled 
by corresponding contributions contained in N^.^{E). This is demonstrated in pT| . 

Finally, the spectral determinant can be factorized in the form A{E) = A(0) Di^t{E) D^xtiE) 
where Di^t{E) and D^xtiE) are completely specified in terms of the inside and outside billiard 
problem, respectively |28|] . 

In the present paper we do not aim at a generalization of the two-dimensional results to 
three dimensions. We rather use the boundary integral equation as starting point for the 
derivation of semiclassical contributions of periodic orbits. One difference in comparison to the 
two-dimensional case is that the integral equation in (|1|) is singular since the integral kernel has 
a l/\f — f'\ singularity. In the surface integral (|l]) this singularity is integrable, but the trace 
of the integral operator Q is infinite, and the Fredholm determinant has to be regularized. 

This can be done by considering the parameter dependent spectral determinant 

( °° \n ] 
Aa(^) :=det(l-AQ(E)) =exp|-^ — Trg"(E)| , |A| < i? , (8) 

whose representation by the sum over n has a positive radius of convergence R if all traces are 
finite. If the first / traces Q"'{E), n = 1, . . . , /, are infinite one can define a modified Fredholm 
determinant by replacing the sum in (^) by a similar expression in which the first / traces are 
missing ||3^. This expression can be analytically continued to the whole complex A-plane. The 



new determinant again has zeros at the eigenvalues of the Schrodinger equation. In case of the 
integral operator (Rl) the first two traces TtQ and TtQ^ are divergent. 



3 The semiclassical Green function 

We take the derivative of eq. (^ with respect to k, insert (^ for A = 1, and obtain the following 
equation for the level density in terms of the wavenumber k 

1 r/ °° 1 

d{k) = d„{k) + -TTrlm V -TrQ"(A;) , (9) 
vr ak ^ — ^ n 
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where 

TrQ"(A;) = (-2)" / da, . . . da^ df,,Go{f2,n, E) df,,Go{f3,r2, E) . . . df,„Go{n,f,,, E) . (10) 

JdV 

Here the vectors are located on the boundary, dai is an area element at r^, and is the 
outward normal vector at fi. 

As noted in the previous section, all traces of powers of Q in the sum in are finite if the 
sum starts with uq > 3. For the semiclassical approximation this has the disadvantage that 
the oscillatory contributions of n-bounce orbits with n < hq are shifted into the smooth part 
'^sm,no(^) through eq. (j^) with the redefined spectral determinant. This problem can be avoided 
by replacing the integral kernel in @ by its leading semiclassical term. Then the integral kernel 
is non-singular and yields the same leading semiclassical contributions to ([T0| ) as before. For 
the derivation of the semiclassical approximation, i. e. after replacing the integral kernel by its 
leading semiclassical term, no thus can be set to no = 1. 

For real values of /c a second divergence occurs since the sum in is in general divergent 
for real k. One can avoid this by starting with an exact formula for a smoothed level density 



27| . We assume for the following derivation that k has a sufficiently large imaginary part in 
order to make the sum convergent, and that Re A; > 0. 

The integrals that appear in eq. ([lOD can be given a simple interpretation in terms of the 
Green function of the billiard system. According to the multiple reflection expansion of Balian 
and Bloch this Green function can be represented by a sum over partial Green functions 



Q{n) corresponding to n reflections on the billiard boundary, with < n < oo 

oo 

G(f, f ', E)=Y, G(")(r, f E) , (11) 



n=0 

where 



G(")(f, r\ E) = (-2)" / da,... dan Go(ri, r ', E) dn^G^^r^, n, E) . . . dn^r, r;, E) . (12) 

JdV 

A comparison of (p!0[) and (|l^) shows that 

d{k) = 4m,no(fc) - - V - / da dn,G^''-^\r,r' ,E)\,=,, . (13) 

^ dk n Jev 

In the following we derive semiclassical contributions of periodic orbits to d{k) by evaluating 



the integrals in ([T2|) and (0) in stationary phase approximation. The semiclassical evaluation 
of (0) expresses the partial Green function G^""^ in terms of classical trajectories from r ' to r. 
These trajectories can be of three kinds: interior orbits which run inside the billiard domain 
V, exterior orbits which run outside V, and ghost orbits which cross the boundary and run 
inside and outside. However, in the sum ([TT| ) the leading order contributions of ghost orbits 



from different n cancel W^. In the following we restrict to the consideration of interior orbits, 
i. e. we denote by Gic^ the semiclassical contributions to G^"-* from interior orbits. We will show 
that this contribution is given by 

Gi^)(r, f', E) = y" , ^ exp \ikL - i^-v^^ - m\ ^ (14) 

" ^ ^ ^ 47rv/|A;2det5^J ^ I ^" 2 ^" i ' ^ ' 
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where 7„ labels all trajectories that run from r ' to r and are reflected n times on the boundary 
in between. Furthermore, l^^ denotes the length of the trajectory 7„, and u^^ is the number of 
conjugate points from r' to r plus twice the number n of reflections on the boundary. B^^^ is a 
submatrix of the stability matrix for the trajectory. The stability matrix is defined by 

Here dgl ' and dp± ' are two-dimensional vectors which describe infinitesimal, perpendicular 
deviations from a considered trajectory at the starting point. dq± and dp± are the corresponding 
deviations at the end point that are obtained by a linearization of the motion in the vicinity of 
the considered trajectory. Correspondingly, M is a 4 x 4-matrix, and A, B, C and D are 2x2- 
matrices. The determinant of the submatrix B is thus given by det i? = M13M24 — M14M23. 
The denominator in (p^ ) is independent of k since detB is proportional to 1/E. Properties of 
the symplectic stability matrix are given in appendix 0, and a method for its determination is 
described in appendix 0. 

We prove now eq. (|T^) by induction. For n = it is correct since 

Go(f,r-",E)=Gi°)(r,r-",E). (16) 

This follows from the form of Go in (|^) and the fact that z/ = and deti? = (r — r')^ for the 
free motion (see the stability matrix for the free motion in appendix 
It remains to show that 

(-2) / da Gi:^ (f, r-;, E) dnGoin, r, E) ^ G^^^) (f,, r;, E) , (17) 

JdV 

where the approximate sign denotes that the integral is evaluated in stationary phase approx- 
imation. We use the following notation at the point r of the boundary: The momentum of 
the particle for the incoming trajectory is denoted by p', and the momentum of the outgoing 
trajectory is p. Both have modulus k = y/E. 

The normal derivative of the Green function is given in leading semiclassical order by 

dnGoin, f,E)^-ih-p G(°) (ffe, r,E) = zk cos a G^") (n, r, E) , (18) 

where a is the angle between — n and p. The left-hand side of (P!7| ) thus is a sum over terms of 
the form 

f exp{2A;[/W(f,r-;) + /W(n,r)]-^fz/W} 

—Ztk cos a / dsids2 — (19) 

JdV 167r2/(0)^|P det5(«)| 

where the element da has been replaced by dsi ds2. 

The stationary points of the above integral are determined by the conditions 

= ^[^^°Hn,r) + /(")(f,r-;)] = k,2 ■ [-| + j] , (20) 

where ti and t2 are two orthogonal unit tangential vectors at the considered point of reflection. 
In eq. (|19|) and in the following the length of a trajectory is given two arguments when it is 
necessary to specify the start and the end point of a trajectory. The second equation in ( pil| ) 
follows from the relation of the length function to the action S{r,r\E) = kl{f,f'). 
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We neglect the possible solution p = p'of (^) since it leads to contributions of ghost orbits. 



The other solution of (pOD yields the condition for elastic reflection: p, p' and fi lie in one plane, 
the reflection plane, and the angles between —p and n and between n and p' are identical. The 
value of this angle a lies in the range < a; < 7r/2. The sum over all stationary points thus 
expresses the integral J by a sum over all trajectories with n + 1 reflections on the boundary. 

For the following evaluations it will be useful to express the tangential vectors and the 
normal vector in local coordinate systems of the trajectory, in which the first coordinate is 
always in direction of the trajectory, and the other two coordinates are perpendicular to it. 



n = — cos a Ci — sm a 62 
ii = sin a ci — cos a 62 
k = eg 



■ cos a e[ — sin a 62 = cos a e'{ — sin a [cos 6 e'2 + sin 6 63] 

■ sin a e\ + cos a e'2 = sin a e'/ + cos a [cos 9 e'2 + sin 9 63] 



cos 9 63 



sin 6* 62 , 



(21) 



where 12 = n x ii. In (0) the unit vectors e" belong to the local coordinate system for the 
incoming trajectory, and the unit vector e'2 lies in the reflection plane of the previous reflection. 
This coordinate system is rotated around the trajectory such that the new vector e'2 lies in the 
new reflection plane. The unprimed coordinate system is the local coordinate system directly 
after the reflection, and 62 again lies in the reflection plane. 

We continue now with the evaluation of the integral in eq. ( [TPD by stationary phase ap- 
proximation. For that purpose the second derivatives of the exponent have to be determined. 
According to ( PD| ) they consist of two parts, the derivatives of the tangential vectors and the 
derivatives of the momenta. Those of the tangential vectors follow from eq. (|60|) of appendix 
and are given by 



(22) 



dti n dti h dt2 n dt2 fi 

dsi Ra ' ds2 Rc ' dsi Rc ' ds2 Rb ' 

where Ra, Rb are the radii of curvature of the surface in the reflection plane and perpendicular 
to it, respectively. Ra, Rb and Rc are given in ( ^0] ) in terms of the two main radii of curvature 
Ri and R2 and the angle f5 between ti and the main curvature direction for Ri. 

The derivatives of the momenta of the incoming and outgoing trajectory are expressed in 
terms of the local coordinate systems before and after the reflection, respectively. 



dp' _ . 



V)p 
V)p' 



cos a 



cos a 



dp2. , dps^ 
dq2 dq2 
dp'2^, , ^Pa 



dp dp2^ dps^ 

-w- = {t2- V)p = —62 + -Tz—es 
ds2 dqs dqs 



dq'2 dq'2 



3 -/ 



dq's 



dq's 



-3 ' 



(23) 



and with (^), (0) and (p3D the second derivatives of the length function follow as 



d^ 



ai 



hi 



Cl 



di :-- 



dsidsi 

d^ 

ds2dsi 

d^ 

dsids2 

d^ 

ds2ds2 



[/(°)(r-;,r-)+/(")(f,r-;)] 
[/(°)(r-;,r-)+/(")(f,r-;)] 
[Z(°)(r-5„r)+/(")(f,r-;)] 
[/(°)(r-;,r-)+/(")(f,r-;)] 



2 cos a 


+ 


dp2 


dp'2 


cos^ a 


Ra 


dq2 


dq'2_ 


k 


2 cos a 


+ 


'dp2 ^ 
.dq3 


dp'2 


COS a 


Rc 


dq'i. 




k ' 


2 cos a 


+ 


'dp3 ^ 
_dq2 


dp's 


cos a 


Rc 


dq'2. 




k ' 


2 cos a 


+ 


dp3 


dq's. 


1 


Rb 


dqs 


k' 



(24) 
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We discuss in more detail how the partial derivatives of the momenta in p^ ) have to be 
interpreted. The derivatives dpi/dqj describe the change of the component Pi upon a change 
of qj when the other components of q and qh are hold fixed. Thus it can be recognized as 
an element of the A'^-matrix of appendix for the free motion from q to qh. By the relations 
of appendix ^ it can be expressed in terms of the stability matrix Mt for this part of the 
trajectory. Note, however, that primed and unprimed quantities have a different meaning here 
and in appendix^. Here the unprimed quantities are the quantities directly after the reflection, 
i. e. they describe the beginning of the trajectory from q to gj,. In the appendix 0, however, 
the unprimed quantities are those at the end of a trajectory, so the meaning is reversed. 

The derivatives dp[/dq'j on the other hand describe the change of the component p[ upon 
a change of q'j when the other components of q' and are hold fixed. Again these partial 
derivatives are elements of an A^-matrix, but now for the trajectory from to g', including 
the final rotation that transforms the double primed coordinate system into the primed one. 
Again by the relations of appendix ^ they can be expressed in terms of the stability matrix 
M' = Ms M^^\ where the rotation matrix Ms is given in (^81). 

The quantities ai, . . . ,di in can also be expressed in terms of the stability matrix for 
the full trajectory from % to that is given by M := M*^""^^) = Mt Mji Ms M*^"^ where Mji is 
the stability matrix for a reflection given in (|49|). The following expressions can be proved by 
a direct although lengthy evaluation that was done with Maple. One finds that 

[Mi3 - Mi4 M^g] cos^ a [M13 M[^ - Mu M^g] cos a 

" ~ K3 M^, - M[, M^3] ' [M[,M^,-Mi,M^,]l(^' 

[M23 M^4 - M24 My cos a [M23 - M24 My 1 

[M(3 M^4 - M{4 My m ' ' [M[,MU-MUM^,]m- ^ > 

and 

_ [Mi3 M24 - Mi4 M23] cos^g _ cos^ a det B 

' ' ~ [M[, M^4 - M[, M^3] {my ~ {my det b' ' ^ ^ 

where M = Mt Mr M' and M' = Ms M^") as before. With (|2|) the contribution of a stationary 
point to the integral in (0) can be evaluated and results in 

, exp{iA:[/(") + /(o)] + f[ais2 + (6i + ci)siS2 + rfis2]-ifz/(")} 

-'sD = —ik cos a / dsi ds2 = 

J 87r2/(0)^|P det B'\ 

exp \ikm+^^ - - ill] 
= ^ ^ , (27) 

where det-B' = detS*^""^ has been used, and 

if ai > and {ai di — 61 Ci) > , 
^(n+l) ^ ^(n) + 2 + <j 1 if (ai di - 61 Ci) < , (28) 

2 if ai < and (ai di — biCi) > . 



The expression ( pTl) has the correct form of a contribution to the semiclassical Green function 
in ([T^). It remains to show that the condition (|2^) gives the correct index z/^"^^'' that is defined 



as the number of conjugate points from to ri, plus twice the number of reflections in between. 



8 



The index z/*^""*"^) increases in comparison to z/'^") by 2 because of the additional reflection, 
and possibly further by 0, 1 or 2, depending on the number of conjugate points between r and 
fb- To determine these conjugate points we consider deti? as a function of the length l^^^ as 
we vary the end point of the trajectory from r to r?,. A conjugate point occurs at every zero 
of deti?, and there can be at most two zeros between r to Vb, because det-B is a quadratic 
function of This follows from the form of the matrix for free motion Mt in (|^). We note 
that a multiplication of a stability matrix by a reflection matrix changes only the sign of 
detB. 

From these considerations follows that there is exactly one zero of det B between r and Vb if 
(ai di — bi ci) < because then det B has changed sign between r and Vb- If (^i di — hi ci) > 
there can be no or two zeros of detS between r and Vb. This depends on the sign of oi. The 
quantity ai has at most one zero as a function of as can be seen from its dependence on l^^^ 
in 

COS^ Q/ \ 
^1 = ^(0) + ^1 ' ^1 = 62 , Ci = 62 , ^1 = y(5y + ^3 , (29) 

where Ci, 62 and 63 are independent of Furthermore, the dependence of the coefficient 
Oi, . . . , c?! on /"^^^ shows that a zero of Oi is always after a first zero of det S, and before a second 
zero of det B there is always a zero of Oi. (If ai = then (ai rfi — 61 Ci) < and for that reason 
there must be a zero of (oi di — 61 Ci) before since both, ai and (oi di — 61 Ci), are positive for 
small /'•^^ Furthermore, ai and di are monotonously decreasing functions of Z*-*^-*. If aidi = e\ 
for two different values of both, ai and c/i, must have changed sign in between.) Thus the 
number of conjugate points is given correctly by (pSf ) and the semiclassical expression for G*^*^^ 
in (|14|) is correct. 



4 Contributions of periodic orbits 



We continue with the evaluation of the integral (|T^) for the spectral density. Inserting the semi- 
classical expressions (p!8D and (|1^) for the partial Green function into ([T3|) yields contributions 
to the level density of the form 

2 d^ If , exp{zA;/("-i) -zfz^("-i)} 

— — Im — / dsids2ikcosa — , — -. (30) 

-K dk n Jqj, 47r^|A;2 det5("-i)| 



We consider now the contributions of stationary points to the integral in (pOl). The stationary 
points are determined by 

which is the condition for a specular refiection. Again p' is the momentum directly before the 
refiection, and p the momentum directly after the refiection. The contributions thus come from 
periodic orbits with n specular refiections on the billiard wall. We consider one such orbit and 
label it by 7. It has n/r^ distinct refiections on the boundary if is its repetition number. 
Consequently, one obtains n/r-^ identical semiclassical contributions for this orbit, since n/r^ 
different stationary points correspond to it. (That the contributions are identical follows from 
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the fact that a stationary phase contribution to Tr in ( [T0| ) does not depend on the order in 
which the n integrals are evaluated.) 

For a determination of the second derivatives of the exponent in ( pOD the derivatives of the 
momenta p and p' are needed. 



dsi 



ds-2 



dsi 



dsi 





V)p 


ik 


V)p 




V)p' 




V)p' 


ts in (121) 



dp2 . dp3 dp2 . dp3^ 
Oq2 oq2 oq, oq, 

dp2. dp3^ dp2^ , dp3^ 

dgs dgs dg^ dq'^ 



cos a , 



5g2 



dq- 



-e, + 



5^2 ^ ^^2 



COS a , 



r, -62 + Tf'eg + -^62 + -^e;, . (32) 
dq^ dqs dq's dq'^ 

^ there are additional terms in ( |32D . They arise because a change 

of the reflection point changes the initial point as well as the end point of the trajectory from 
fto f. With (|31|), dH) and (H) we obtain the following expressions 



a2 



C2 



do 



ds2dsi 
dsids2 
dsi 



[/{n-l)| 
[/{n~l)| 



r, r 



r, r 



r, r 



r, r 



)] 



2 cos a 


+ 




5^2 


_ dp'2 




cos^ a 


Ra 




dq2 




k 


2 cos a 


+ 


.<9g3 


9p2 ^ 

dq's 


M + 
dqs 


dp2 


cos a 




dq's. 




k ' 


2 cos a 


+ 


'dp3 _ 


dps _ 


dPs 1 

5g2 


dPs' 


cos a 


Rc 


_dq2 


dq'2 


dq',_ 




k ' 


2 cos a 


+ 


dps 


dps 


5g3 


dq's. 


1 


Rb 


dq3 


dq's 


k' 



(33) 



Again all partial derivatives of the momenta are elements of an A^-matrix that corresponds to 
the stability matrix M' = Ms M^^~^\ and thus the quantities 02, ... ,^2 can be expressed in 
terms of the elements of M'. They can further be expressed in terms of the full stability matrix 
of the orbit M = = M^i M', where describes the reflection. One finds that 



02 



C2 



M 



42 



kdetB 

Ms2 

kdetB 



cos^ a 



cos a 



M 



41 



kdetB 

Msi 
kdetB 



cos a 



and 



a2 d2 — &2 C2 



det(M - l)cos2a 



(34) 



(35) 



The quantity A4ij in (|3 



fc2 det B 

is the determinant of the 3 x 3-matrix that is obtained by dropping 
the 2-th row and the j-th column of the matrix (M — 1). 

With these expressions a stationary phase approximation for the contribution of an isolated 
periodic orbit results in 

/("-i)cosa. /•, , fcexp{2fc/("-i) + f(a2S? + (62 + C2)siS2 + o?2si; 



d^{k) 



Im / dsi ds-2 



^y\k^ det5("-i)| 



vrr^ ^|det(M^ - 1)| 



Reexp 



(36) 
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where the relations I det 5*^" ^^1 



det B\, 



/("-I) and 



if 02 > and (02 d2 — 62 C2) > , 
/i-y = //^""'^ + 2 + <{ 1 if (a2 c/2 - &2 C2) < , 

2 if 02 < and (02 (^2 — 62 C2) > , 



(37) 



have been used. Eq. (^) is the expected contribution of an isolated periodic orbit. The relations 
for determining the Maslov index /i of a periodic orbit are summarized in appendix y. The 
conditions for the additional contributions to /i^ in ( pT]) can be shown to be identical to the 



general conditions in . 



In applications three-dimensional billiard systems often have an axial symmetry. Then 
periodic orbits appear in one-parameter families. In the following we derive the contribution of 
such a family. A general treatment of semiclassical trace formulas in the presence of continuous 
symmetries can be found in [^. In this article, a general formula is given for families of orbits 
in axially symmetric systems. The formula which is obtained in the following has a different 
form which is not explicitly independent of the starting point of the trajectory, but it has the 
advantage that it is expressed directly in terms of the stability matrix M. Of course, both 
formulas must yield the same result. 

In presence of an axial symmetry the quantity (02 ^2 — &2 C2) vanishes and the symmetric 
matrix {"^2%) eigenvalues and (02 + ^2). By a rotation of the coordinate system the 
matrix can be diagonalized and one obtains 



/("-I) cosa , /•,,,, fcexp + f (a2 + d2)s'i - ^'^'/"-i) 

Im / — 



27r2r, 



v/|A;2 det5("-i)| 



l-fp cos a 



Re 



27Tk 



, , . . o , -n — 7 exp < ikL 

\kM42Cos^a + kM3i\ 1 ^ 2^^ 4 



(38) 



where 1^ = 1^'^ and 



(n-1) 



+ 2 + 



if (a2 + c/2) > , 

1 if (02 + c/2) < . 



(39) 



In (|38|) the integral over S2 was evaluated over the range 211 p/N where p is the distance 
of the stationary point on the billiard surface to the symmetry axis, and is the number of 
distinct rotations around the symmetry axis that map the orbit onto itself. The result (^) does 
not depend on the reflection point that is chosen as starting point of the trajectory, although 
this is not explicit from its form. The reason is again, that a semiclassical contribution to Tr Q"^ 
in (|10D does not depend on the order in which the n integrals are evaluated. 



We acknowledge financial support by the Deutsche Forschungsgemeinschaft under contract 
No. DFG-Ste 241/6-1 and /7-2. 
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A General properties of the stability matrix 



The stability matrix M is defined by the following relation 



/dg2\ 

\dp3J 



M 



dp'2 



M 



/Mil M12 Mi3 MiA 

M21 M22 M23 M24 

M31 M32 M33 M34 

\M4i M42 M43 M44/ 



(40) 



Here dgg, dgg, dp2, and dpg are infinitesimal deviations from the considered trajectory at the 
starting point. They are given in a local coordinate system in which the first coordinate is along 
the trajectory, and the other two coordinates are perpendicular to it. The unprimed quantities 
dq2, dqs, dp2 and dp^ are the corresponding infinitesimal deviations at the end point of the 
trajectory. They follow from a linearization of the motion in the vicinity of the considered 
trajectory. In a compact notation 



dq_i 
dpA 



M 



dq_i_ 
dp± 



M 



A 
C 



B 
D 



(41) 



where A, B, C and D are 2 x 2- matrices. 

The stability matrix is a symplectic matrix. It satisfies M'^ J M = J, where J denotes the 
matrix J = ( i*/ ) ^ is the two-dimensional unity matrix. The inverse of M is given by 
= —JM^ J, and from M = M = 1 one obtains the following relations between 
the submatrices of M 



A^C 
AB^ 



C^A 
BA"^ 



B^D 
CD^ 



D^B 
DC^ 



A^D 
AD^ 



C^B 
BC^ 



I . 



(42) 



Both rows of equations are equivalent to each other and to the symplectic condition for M. 
Altogether they yield 6 independent conditions for the matrix elements of M, so that only 10 
elements of M are independent. 

Another useful quantity is the matrix which is defined in the following. It is obtained by 



solving the linear equations in (|41|) for pj_ and p±_ ' 



dpx' 
dpi 



A^ 



dgl' 
dgl 



A B\ /dgl' 
C DJ Ugl 



(43) 



The relations between the submatrices of A^ and those of M are given by: A = —B^^A, 
B = B^^, C = C — D B^^A, and D = D B^^. On the other hand, the matrix elements of A^ 
can be obtained from the action function of the trajectory S{q, q', E) through 



AT, 



cr,- 



dzidzj 



dgl' 
dgl 



-1 



1,2 
3,4 



(44) 



Since the mixed second derivatives of the action function do not depend on the order of the 
derivation this yields 6 conditions for the matrix elements of A^ which are equivalent to A = A^, 
B = -C^, and D = D^. 
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B Determination of the stability matrix 



In two-dimensional billiard systems the stability matrix is composed of two kinds of partial 
stability matrices Mt and where Mt denotes the matrix for the free motion between two 
reflections and is the matrix for a reflection. They have the following form (see [^]) 



L/p 
1 



Mr 



-1 

_Rcos a 







(45) 



where L is the length of a trajectory between two reflections, R is the radius of curvature at a 
reflection point, a is the angle of incidence, and p is the modulus of the momentum. 

In three dimensions there is an additional kind of matrix which describes a rotation of the 
local coordinate system around the trajectory The stability matrix for a trajectory from 
a point fa to a point rf, with n reflections on the boundary can be written in the following form 



M 



'Ml M^ . 



M^r 



■2 Ml Ml 



Ml Ml M]: 



(46) 



Here Mt is the matrix for a part of the trajectory of length L between two reflections. It 
is a straightforward generalization of the matrix Mt in two dimensions. Ms corresponds to a 
rotation of the local coordinate system around the trajectory such that the new coordinate with 
index 2 lies in the reflection plane that is spanned by the incoming and outgoing trajectory at 
a reflection point. The angle of rotation 9 is defined by 

62 = €-2 COS 6^ + 63 sin 6* , dg2 = dgg cos 6 + dgg sin 6 , 

63 = —62 sin 6* + 63 cos , dg3 = —dq'2 sin 6 + dg3 cos 6 , (47) 

where e[ are the unit vectors along the old coordinate axes and Cj are the unit vectors along 
the new coordinate axes. Mt and Ms are given by 

/ cos 6 



Ms 



(I L/p \ 

1 L/p 

1 

\0 1 / 

Furthermore, Mji is the matrix for a reflection 

/ -1 


Mn = 2p 



sin 6 
sin 6 cos 6 





Ra COS a 
2p 

Rr 




1 

2p 

Rc 
2p cos a 

Rb 
















o\ 








-1 








V 



\ 



COS 9 sin 9 

— sin 9 cos 9 J 



(48) 



(49) 



Here a is the angle of incidence, i. e. the angle between the incoming (or outgoing) trajectory 
and the normal to the boundary. Ra and Ri, are the radii of curvature in the reflection plane 
and perpendicular to it, respectively. Ra, Rb and Rc can be expressed in terms of the two main 
radii of curvature at the reflection point -Ri and R2, and the angle /3 between the tangent lying 
in the reflection plane and the direction of the main curvature 1/Ri (see eq. (|57D). 

1 _ cos2 p ^ sin^ p 1 _ sin^ /? ^ cos^ /? 1 _ R2 - Ri 

Ra Ri R2 ' Rb Ri R2 ' Rc R1R2 



COS /3 sin /? . (50) 
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The form of the matrices M^, Ms and Mr has been given in [0 for the case of reflections 
on spheres [Ri = R2). In ||3^ these results were generalized to cases where one of the main 
curvature directions of the billiard boundary lies in the reflection plane {6 = 0). Here we give 
the form of Mr for arbitrary reflections. The derivation is carried out in appendix 0. 

Finally, we discuss the stability matrix for a periodic orbit. It is given by ( ^BD if a is replaced 
by n and the term is dropped (see eq. (0)). In this form it is assumed that the local 

coordinate system after one traversal of the trajectory has the same orientation as it has at the 
beginning. Otherwise an additional rotation matrix Ms has to be multiplied to the product. 
However, since one has the freedom to choose the orientation of the local coordinate system at 
the starting point, one can always choose it to have the same orientation as at the end point. 
In other words, for the evaluation of (|5lD from right to left one starts with a local coordinate 
systems whose second coordinate lies in the reflection plane at fVi- 



C Determination of the Maslov index 

We summarize in this appendix the results of sections ^ and ^ for the Maslov index of a 
trajectory. Consider a trajectory with n reflections on the billiard wall. If one of the reflection 
points is chosen as starting point of the trajectory the stability matrix can be written in the 
following form 

M = M^M^ M^^"-i . . . M|^2 j\^2 ^2 ^2^1 ^1 ^i^n _ (-5^^ 

Here the orientation of the local coordinate system at the initial point has to agree with its 
orientation at the final point as discussed in the appendix 0. Now the trajectory is traversed 
once, and correspondingly the multiplication of partial matrices in (^) is carried out. The 
Maslov index changes by two for every reflection matrix in this product (in case of Dirichlet 
boundary conditions). Otherwise it changes only during a part of the trajectory between two 
reflections. How much it changes is determined by the following method. For every M^^-matrix 
in (^) except for the first one M^*~" one denotes the product of partial matrices up to this 
Mr-matrix by M', and the product including this M^-matrix by M such that M = Mt M' . 
Then one considers the following two quantities 

_ Mi3 M^4 - Mi4 M^3 _ M13M24-M14M23 . . 

' M(3M^4-M(4M^3' ' M{3M^4-M{4M^3- ^ > 

The Maslov index changes according to the following rule 

{0 if Zi > and Z2 > , 
1 if < , (53) 
2 if Zi < and Z2 > , 

where v' denotes the previous value of the index. This follows from (^) and (p8|) by noting 
that the Mij-element change only their sign upon a reflection and the M2relements remain 
unchanged. (Note that M' is defined here including the reflection matrix in contrast to the 
main section.) 

Finally there can be an additional change of the Maslov index that arises from the evaluation 
of the trace. It is determined by considering 

7 -^42 „ det(M-l) 
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where M is now the total stabihty matrix of the periodic orbit given in (0), and B is its upper 
right submatrix as defined in appendix 0. A^^ is the determinant of the 3 x 3-matrix that is 
obtained by dropping the i-th row and the j-th column of the matrix (M — 1). 
Then the index v is changed according to 

{0 if Zi > and ^2 > , 
1 if ^2 < , (55) 
2 if Zi < and Z2 > . 

This yields the total Maslov index /i of the periodic orbit. 



D Derivation of the stability matrix for a reflection 

Let the point of reflection be at the origin of the coordinate system. The x'-^- and of the 

coordinate system are chosen to lie in the directions of the two main curvatures of the boundary 
surface. Then the boundary is locally represented by 

/ 2 / 2 

where Ri and R2 are the two main radii of curvature. Now the coordinate system is rotated 
about the Xg-axis so that the new Xi-axis lies in the plane of reflection, and the component of 
the incoming momentum in xi-direction is non-negative 

x[ = xi cos /5 + X2 sin /3 , x'2 = —xi sin P + X2 cos P , x'^ = x^ . (57) 

In the new coordinates the boundary is locally represented by 

2 2 

_X]_ X2 X1X2 , . 

ZKa ZKb Kc 

where Ra and Rb are the radii of curvature in direction of the new coordinate axes in planes 
containing the normal vector, and the relations between Ra, Rb, Rc and Ri, R2 and /3 are given 
in m. 



Quantities before the reflection are given a prime in the following, and the conserved absolute 
value of the momentum is denoted by p. The angle of incidence a is defined as the angle 
between the incoming momentum p' and the outward pointing normal vector to the boundary 
n, < a < 7i/2. The outgoing momentum is p. One has 



p = p' — 2{p' ■ n)h = p \ . (59) 





0^ 




1 sina \ 


n = 


V ? } 


, p =p\ 











V cos a / 




In the following one considers an infinitesimal change in the position and the momentum of the 
incoming trajectory and determines the corresponding deviations for the outgoing trajectory. 
First we discuss how the normal vector changes if the point of reflection is changed. The new 



15 



normal vector is determined from two tangential vectors at the point r = (xi, 0:2, ^3) with X3 

( \ 



given by eq. (^) 



df 



1 




Xi X2 
\ Ra RcJ 



df 

8X2 




1 



X2 



Xi 

RcJ 



(60) 



from which one obtains 



n 



h X t2 
\ti X fsl 



Xi 


+ 


X2 


Ra 


Rc 


X2 


+ 


Xi 


Rb 


Rc 



( 



\ 1 / 



Ra Rc 



X2 _^ 

Rb Rc 



1 -1/2 



(61) 



Now one considers the four cases that one of the infinitesimal quantities dq'2, dgg, dp^ and dpg 
is different from zero. They describe infinitesimal deviations from the incoming trajectory in 
position and momentum, respectively, immediately before the reflection. They are expressed in 
a local coordinate system, where the first coordinate is parallel to the trajectory and the other 
two coordinates are perpendicular to it. The quantities with index 2 are in the plane of reflection 
and those with index 3 perpendicular to it. For each of these four cases the corresponding 
deviations for the outgoing trajectory immediately after the reflection are determined. They 
are written without a prime. 



dp2 7^ 0: One sees easily that only dp2 is different from zero: dg2 
and dp3 = 0. 



0, dga = 0, dp2 



-dp'2 



dpg 7^ 0: Now only dpa is different from zero: dg2 = 0, dgs = 0, dp2 = and dpa = dpg, 



dg2 7^ 0: The deviation of the point of reflection is given by dxi = dgg/ cos a = — dg2/ cos a 
and dx2 = 0. It follows that dg2 = — dgg and dgs = 0. With the normal vector in (|6TD , 
and with p' and the relation for p in (K^) one finds: 



n 







Ra COS a 






, p = p 


Rc cos a 




\ 1 > 





( 



sma 



2dg^ \ 



R„ 



\ 



cos a. — 



2dq^ 

Rc 

2 tan a dgg 

Rn 



( 



Xdp2 \ 



7 



p sm a — cos a 
dps 

y— p cos a — sin adp2 j 



(62) 



From this follows that dp2 = 2pdg2/(-Ra cosa) and dps = —2pdq2/Rc- 

dgg 7^ 0: The deviation of the point of reflection is given by dxi = and dx2 = dgg = dgg, and 
it follows that dg2 = and dgg = dg3. With the normal vector in (pl|), and with p' and 
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the relation for p in (p9[) one finds 



n 



Rb 

V 1 / 



p = p 



/ . 2 cos a dq'n \ 
' sm a > 

Rc 

2 cos a dgg 



V" 



COS a 



Rb 

2 sin a dgg 

Rr 



7 



/ p sin a — cos a dp2 \ 
y — p cos a — sin a dp2 y 



From this follows that dp2 = '^pdqyRc and dps = — 2pcosa dgg/i?;,. 
Collecting all results one obtains the form of the stability matrix for a reflection in 



(63) 
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